#################################################################
#bvn function set up the expression of the bivariate normal dstn
#function
#input:              x - vector x
#                    y - vector y
#                    mx- mean(x)
#                    my- mean(y)
#                    sx- stdev(x)
#                    sy- stdev(y)
#                    rho-correlation(x,y)
#output              evaluated BVN function value for given data (x,y)
#################################################################

bvn<-function(x,y,mx,my,sx,sy,rho) {

#form the var-cov matrix
sigma<-matrix(c((sx*sx),(rho*sx*sy),(rho*sx*sy),(sy*sy)),nrow=2)
#form the matrix of X-mu_X
xv<-matrix(c(x-mx,y-my),nrow=2)
#evaluate the function
eval_f=1/2/det(sigma)^.5*exp(-.5*t(xv)%*%solve(sigma)%*%xv) 
#return evaluatet function value
return(eval_f)

}
#################################################################


#################################################################
#bvngrid function set up the grid for ploting the values obtained 
#from bvn function
#input:              x - vector x
#                    y - vector y
#                    mx- mean(x)
#                    my- mean(y)
#                    sx- stdev(x)
#                    sy- stdev(y)
#                    rho-correlation(x,y)
#output              matrix containing the values of evaluated bvn 
#                    function over the grid

#################################################################

bvngrid<-function(x,y,mx,my,sx,sy,rho) {
  
#initialize empty matrix z with dimension
# length(x) x length(y)
z <- matrix(NA,length(x),length(y))
#fill the matrix z with evaluated function values
for (i in 1:length(x)) {
for (j in 1:length(y)) {
z[i,j] <- bvn(x[i],y[j],mx,my,sx,sy,rho) }}
#retun matrix z
return(z)
}
#################################################################



#mesh plot and contour plot usualy go together 
#seeing both of them give us better understanding of
#the geometry of the distribution


#################################################################
#set up the values of (x,y) for the grid
x<-y<-seq(-3,3,by=0.1)
#evaluate matrix z
z<-bvngrid(x,y,0,0,1,1,.5)
#draw a mesh (function surface) plot
persp(x,y,z,theta=25,phi=-5)
#make a contour plot
contour(x,y,z)
#################################################################

#theta=20 and phi=-5 parameters of the persp plot
#are for seeing the mesh plot from different angles
#play with these values in order to see the function surface
#from different angles

#################################################################
#draw another plot with different mu and Sigma
#set up the values of (x,y) for the grid
x<-y<-seq(-3,3,by=0.1)
#evaluate matrix z
z<-bvngrid(x,y,0,0,1,1,0)
#draw a mesh (function surface) plot
persp(x,y,z,theta=20,phi=-5)
#draw a contour plot
contour(x,y,z)
#################################################################



#################################################################
#draw another plot with different mu and Sigma

#set up the values of (x,y) for the grid
x<-y<-seq(-3,3,by=0.1)
#evaluate matrix z
z<-bvngrid(x,y,0,0,1,1,0.8)
#draw a mesh (function surface) plot
persp(x,y,z,theta=20,phi=-5)
#draw a contour plot
contour(x,y,z)

#################################################################